Three-Dimensional Finite-Element Code for Electerosurgery and Thermal Ablation Simulations

ABSTRACT

A method for performing a simulation of an electrosurgical procedure is provided. The method includes the step of generating a three-dimensional model of a medium comprising one or more material regions representing a physical object or a region of the tissue. The method also includes the step of calculating one of an electrical energy solution and a thermal transport solution based on one of a mathematical model and finite element method. The method further includes the step of updating one or more property values of the material region in response to temperature changes based on an empirically derived value, wherein one of the electrical energy solution and the thermal transport solution is updated as a function of the updated property value.

CROSS-REFERENCE TO A RELATED APPLICATION

This application is a continuation of a U.S. application Ser. No.11/314,332 filed on Dec. 21, 2005 which claims priority to U.S.Provisional Application Ser. No. 60/638,560 filed on Dec. 23, 2004, theentire contents of which are hereby incorporated by reference herein.

BACKGROUND

1. Field

The present disclosure relates generally to data processing and computersimulation systems, and more particularly, to a method and system forsimulations electrosurgical RF thermal ablation processes.

2. Description of the Related Art

Electrosurgery procedures involve the application of RF electric fieldsto produce local heating. The goal is to destroy or to alter selectedtissues. Applications include blood-vessel cauterization, incisionsealing and minimally-invasive procedures for tumor destruction. Bytheir nature, electrosurgery procedures induce large changes in thethermal and electrical properties of tissue. In many procedures, theapplication of RF power follows a complex temporal sequence to ensureeffective local treatment while minimizing damage to surroundingtissues. In the past, the design of equipment and the choice ofoperation sequence relied primarily on empirical experience. As aresult, investigation of new applications is generally a long drawn out,difficult and costly process.

Computer simulations can play an important role in reducing the cost andtime associated with certain experiments. In other words, simulationsoffer a relatively quick and inexpensive approach to fill in missingdata points, to investigate novel methods or to optimize existingprocedures.

Considering the complexity of electrosurgical procedures, genericfinite-element software is generally not very useful. Conventionalelectrosurgical simulation software relies on two-dimensionalsimulations in cylindrical or planar geometries. Such software suitesare incapable of determining end effects for many electrosurgicalprocedures. These software suites are also unable to provide solutionsfor simple electrosurgical probes which may be significantly altered bythe presence of nearby asymmetric structures like blood vessels.

Therefore there exists a need for a system and method for simulation ofelectrosurgical procedures which provides solutions in athree-dimensional medium.

SUMMARY

A method and system for simulating the effect of electrosurgicalprocesses such as RF tumor ablation are provided. The method accordingto the present disclosure covers the complete calculation process frommesh generation to solution analysis. The method employsthree-dimensional conformal meshes to handle cluster probes and otherasymmetric assemblies. The conformal-mesh approach is advantageous forhigh-accuracy surface integrals of net electrode currents. Thesimulation algorithm performs coupled calculations of RF electric fieldsin conductive dielectrics and thermal transport via dynamic solutions ofthe bioheat equation. The boundary-value RF field solution is updatedperiodically to reflect changes in material properties.

The simulation algorithm also features advanced material models with theoption for arbitrary temperature variations of thermal and electricalconductivity, perfusion rate, and other quantities. The algorithmhandles irreversible changes by switching the material reference ofindividual elements at specified transition temperatures. The simulationalgorithm is controlled through a versatile interpreter language toenable complex run sequences. The algorithm can also automaticallymaintain constant current or power, switch to different states inresponse to temperature or impedance information, and adjust parameterson the basis of user-supplied control functions.

According to one embodiment of the present disclosure a method forperforming a simulation of an electrosurgical procedure is disclosed.The method includes the step of generating a three-dimensional model ofa medium including one or more material regions representing a physicalobject or a region of the tissue. The method also includes the step ofcalculating one of an electrical energy solution and a thermal transportsolution based on one of a mathematical model and finite element method.The method further includes the step of updating one or more propertyvalues of the material region in response to temperature changes basedon an empirically derived value, wherein one of the electrical energysolution and the thermal transport solution is updated as a function ofthe updated property value.

According to another embodiment of the present disclosure, a system forperforming a simulation of an electrosurgical procedure is disclosed.The system includes one or more processors coupled to acomputer-readable media configured to store a set of instructionscapable of being executed by the at least one processor. Theinstructions include the steps of generating a three-dimensional modelof a medium including one or more material regions representing aphysical object or a region of the tissue. The instructions also includethe step of calculating one of an electrical energy solution and athermal transport solution based on one of a mathematical model andfinite element method. The instructions further include the step ofupdating one or more property values of the material region in responseto temperature changes based on an empirically derived value, whereinone of the electrical energy solution and the thermal transport solutionis updated as a function of the updated property value.

According to a further embodiment of the present disclosure, acomputer-readable storage medium storing a set of computer-executableinstructions for performing a simulation of an electrosurgical procedureis disclosed. The computer-executable instructions include the step ofgenerating a three-dimensional model of a medium including one or morematerial regions representing a physical object or a region of thetissue. The instructions also include the step of calculating one of anelectrical energy solution and a thermal transport solution based on oneof a mathematical model and finite element method. The instructionsfurther include the step of updating one or more property values of thematerial region in response to temperature changes based on anempirically derived value, wherein one of the electrical energy solutionand the thermal transport solution is updated as a function of theupdated property value.

BRIEF DESCRIPTION OF THE DRAWINGS

The above and other aspects, features, and advantages of the presentdisclosure will become more apparent in light of the following detaileddescription when taken in conjunction with the accompanying drawings inwhich:

FIG. 1 is a schematic diagram of an exemplary computer system;

FIG. 2A shows a non-conformal mesh including box elements;

FIG. 2B shows a conformal mesh including hexahedron elements;

FIG. 2C shows a three-dimensional hexahedron element;

FIG. 3 is an exemplary generated surface of a medium;

FIG. 4 is a graph defining a temperature-dependent electricalconductivity;

FIG. 5 is a screen shot of the software of the present disclosure forradio frequency electric fields illustrating analysis capabilities in atwo-dimensional slice plot mode;

FIG. 6A illustrates an isothermal surface with gray scale codingdetermined by the magnitude of thermal flux;

FIG. 6B is a representation of region boundaries of the probe and sheathgray scale coded by temperature;

FIG. 7 shows voltage variation required to maintain constant currentwith decreasing conductivity for two values of a drive current;

FIG. 8 shows a two-dimensional slice plot of the electrical solution;

FIG. 9A is a plot of probe impedance over time in a solution volume;

FIG. 9B is a plot of maximum temperature over time in the solutionvolume; and

FIG. 10 shows spatial variation of electrical conductivity in aparticular plane; and

FIG. 11 shows a flow chart illustrating a method according to thepresent disclosure.

DETAILED DESCRIPTION

Preferred embodiments of the present disclosure are described hereinwith reference to the accompanying drawings. In the followingdescription, well-known functions or constructions are not described indetail to avoid obscuring the present disclosure in unnecessary detail.

1. Introduction

It is to be understood that the present disclosure may be implemented invarious forms of hardware, software, firmware, special purposeprocessors, or a combination thereof. In one embodiment, the presentdisclosure may be implemented in software as an application programtangibly embodied on a program storage device. The application programmay be uploaded to, and executed by, a computer system 100 including anysuitable architecture such as a personal computer, a workstation orserver. The system 100 is implemented on a computer platform havinghardware such as one or more central processing units (CPU) (e.g.,processor 102), a random access memory (RAM) 104, a read only memory(ROM) 106 and input/output (I/O) interface(s) such as a keyboard 108,cursor control device (e.g., a mouse or joystick) 110 and display device112. A system bus 115 couples the various components and may be any ofseveral types of bus structures including a memory bus or memorycontroller, a peripheral bus, and a local bus using any of a variety ofbus architectures. The computer platform also includes an operatingsystem and micro instruction code, which may be stored in a data storagedevice 114. The various processes and functions described herein mayeither be part of the micro instruction code or part of the applicationprogram (or a combination thereof) which is executed via the operatingsystem. In addition, various other peripheral devices may be connectedto the computer platform by various interfaces and bus structures, suchas a parallel port, serial port or universal serial bus (USB), forexample, additional storage devices and a printer.

It is to be further understood that because some of the constituentsystem components and method steps depicted in the accompanying figuresmay be implemented in software, the actual connections between thesystem components (or the process steps) may differ depending upon themanner in which the present disclosure is programmed. Given theteachings of the present disclosure provided herein, one of ordinaryskill in the related art will be able to contemplate these and similarimplementations or configurations of the present disclosure.

The simulation algorithm according to present disclosure may be run onpersonal computers under Microsoft Windows™ or Linux™ operating systemsand/or software. The algorithm includes an interactive mesh generator,solution programs for RF fields and thermal transport and dedicatedpost-processors for plots and data analysis. The method according topresent disclosure utilizes conformal hexahedron meshes forhigh-accuracy surface integrals, calculates electrical and thermalsolutions with temperature-dependent effects, incorporates advancedmaterial models, and further includes a built-in script interpreter witha flexible control language for autonomous operation.

The following section, Section 2, describes the nature of conformalhexahedron meshes and the tools used for their creation. Section 3covers solution methods for coupled RF electric fields and thermaltransport. Section 4 reviews techniques for modelingtemperature-dependent properties of materials (including phase changes)while Section 5 describes methods to operate probes at fixed current oras part of a temperature-control loop. Features of the interpreterlanguage for sequence control are described in Section 6. Features ofthe postprocessors are covered in Section 7. Section 8 presents examplesthat illustrate some advanced features of the algorithm.

2. Mesh Generation

The foundation of the finite-element method for field solutions is todivide the solution space (e.g., a segment and/or medium 8 of asimulated surgical site) into small volumes (e.g., elements). Thedivision is performed in such a way that each element is uniquely partof one material region. In particular, the material regions may alsorepresent electrodes, insulators or different types of tissue (e.g.,blood vessels, muscles, etc.). This process simplifies the calculationsas well as the organization of the code—only one material model isapplied to each element. The division facilitates the conversion of thegoverning partial differential equations for RF fields and temperatureto large sets of coupled linear equations. Such equation systems areeasily solved on a digital computer, an example of which is discussedabove.

The simplest way to divide space, e.g., medium 8, is the non-conformalmesh 10 of FIG. 2A which employs box elements. However, this approachleads to a poor representation of surfaces between material regionssince curved surfaces are approximated and are not accuratelyrepresented. As a result, field calculations near the surface ofmaterial regions are inaccurate. For example, surface integrals ofcurrent over the sphere represented by non-conformal mesh may beinaccurate.

FIG. 2B shows a conformal mesh 12 representing the medium 8, which giveshighly-accurate results in surface field calculations. In a conformalmesh, elements are flexed near region boundaries so that facets lieclose to the desired shape. The box elements are transformed togeneralized hexahedrons. An exemplary hexahedron element 14 is shown inFIG. 2C.

Mesh generation is accomplished by ordering elements 14 according tocubic logic. However, elements 14 may have arbitrary shapes. The element14 includes a plurality of facets 16 and nodes 18 (e.g., indices). Thefacets 16 may meet at a variety of angles allowing for flexible shapeswhich more closely approximate the medium. In other words, each elementhas a unique set of indices along the coordinate axes such that element(I+1, J, K) is adjacent to element (I, J, K). Mesh structure allowsfaster field solutions, where the boundary-value solution for RFelectric field may be updated hundreds of times during execution withtemperature-dependent materials.

FIG. 3 illustrates an exemplary generated surface of the medium 8, whichincludes mathematically-generated surfaces from a Table 1 containingelevation values. The mesh in the FIG. 3 combines a mathematical surfacerepresentative of tissue 20 with an extrusion (e.g., electrosurgicalprobe 32).

Table 1 shows the script to generate the shapes. In the presentdisclosure, the program was adapted to use constructive solid geometrytechniques to fabricate complex assemblies from simple parts (i.e.,cylinders, spheres and cones, extrusions, turnings, etc.).

TABLE 1 MetaMesh script to generate FIG. 3 PART Type Surf 3DFIG.D02.DATE Name TestSurface Region TestSurface Fab −1.00 SurfaceRegion SVolume 1.00 END PART Type Extrusion L −1.0 −1.0  1.0 −1.0 S L 1.0 −1.0  1.0  0.0 S A  1.0  0.0  0.0  1.0 0.0 0.0 S A  0.0  1.0 −1.0 0.0 0.0 0.0 S L −1.0  0.0 −1.0 −1.0 S End Region Electrode Fab 9.00Surface Region SVolume Rotate 0.0 0.0 180.0 END

3. Field Solution Techniques

Electrosurgery simulations involve two finite-element solutions:radiofrequency (“RF”) electric fields and thermal transport Structuresgenerated in the medium are common to both solutions. Consequently, themethod of the present disclosure uses the same mesh for the RF andthermal solutions, assigning different physical properties to theregions. Because there is little change in thermal properties over aperiod of the RF field, the electric field solution is handled as aboundary value problem with periodic updating. The following equationgoverns the RF electric field solution:

$\begin{matrix}{{\nabla{\cdot \left\lbrack {\left( {ɛ_{r} - \frac{j\sigma}{ɛ_{0}\omega}} \right){\nabla\varphi}} \right\rbrack}} = 0} & (1)\end{matrix}$

In Equation 1, ∈_(r) is the relative dielectric constant of the element14 and σ is the electrical conductivity. The quantity φ (thequasi-static potential) is a complex number. The amplitude and phase ofthe electric field is given by:

E=−∇φ  (2)

The electric field provides one of the source terms for the thermaltransport equation.

q _(e) =σE·E*/2  (3)

Equations 1 and 2 include effects of both real and displacement currentsand hold in the non-radiative limit. The validity condition is that theradiated electromagnetic energy leaving the patient is small compared tothe deposited energy. The condition holds when the RF wavelength islarge compared to the treatment area.

Conversion of Equation 1 to finite-element form uses the standardminimum-residual method with quasi-linear shape functions for thehexahedron elements 14. The resulting linear equation set determines φat the element 14 node locations. Each equation involves coupling to 27neighboring nodes. The coupling coefficients are computed by takingnumerical integrals over elements using the normal coordinate method.The equations are solved with an interactive matrix-inversion technique.An initial electric field solution with 400,000 nodes requires about 2minutes on a 3.0 GHz computer. The solution time is shorter for updateswith temperature-dependent ∈_(r) or σ if changes are relatively small.

The method of the present disclosure solves the following relationship(often known as the bioheat equation) for thermal transport;

$\begin{matrix}{{{pC}_{p}\frac{\partial T}{\partial t}} = {{\nabla{\cdot \left( {k{\nabla T}} \right)}} + q_{e} + q_{m} + {W_{b}\left( {T_{b} - T} \right)}}} & (4)\end{matrix}$

In Equation 4, ρ is the density of the medium, C_(p) is the specificheat, k is the thermal conductivity, q_(e) is the resistive source fromEquation 3 and q_(m) is any other source such as metabolism. The termW_(b)(T_(b)−T) is a thermal source or sink resulting from thevolume-averaged perfusion of blood through capillaries. Here, W_(b) isthe blood flow rate and T_(b) is the incoming blood temperature.

The system and method of the present disclosure recast theminimum-residual formulation for the conformal mesh to reduce couplingto the six nearest nodes. The resulting solutions are fast and stablefor all mesh geometries and choices of time step at the expense of somereduction in accuracy.

Execution of the algorithm consists of an initial RF solution followedby a thermal solution over a specified time span. Temperature-dependentthermal quantities are updated each time step. The electric field mustbe updated periodically if materials have temperature-dependent valuesof ∈_(r) or σ. Electrical field updates may consume considerable time,so it is essential to minimize operations in three-dimensional solutionswith large meshes. Accordingly, the code implementing the method hasbeen carefully analyzed to eliminate unnecessary or redundantoperations. Several techniques have been implemented to reduce run time,such as variable time steps with automatic adjustment are used in thethermal solution. Further, whenever possible, the RF solution is scaledrather than recomputed. The program utilizes precomputed baseline RFsolutions during pulsed RF simulations. In addition, element propertiesand coupling coefficients for the RF and thermal solutions are updatedonly for temperature-dependent materials.

4. Material Models and Temporal Variations

Electrosurgery simulations usually involve materials with significantvariations of properties with temperature. Therefore, the system andmethod of the present disclosure, uses a flexible method of tabularfunctions to represent arbitrary temperature and time variations. Atable is simply a text file where each line contains one value of anindependent and dependent variable (e.g., temperature T and thermalconductivity k(T)).

FIG. 4 shows a typical data set to represent a temperature-dependentelectrical conductivity which is shown as a graph of electricalconductivity over temperature. Electrical conductivity of materials(e.g., tissue) varies with temperature, therefore the method and systemaccording to the present disclosure account for these variations byupdating RF field calculations.

The system and method incorporates a software unit to perform cubicspline or linear interpolations. The unit is designed so that tablesappear as a continuous function to the calling program. The routineshandle non-uniform intervals of the dependent variable. Flexibility isachieved through a free-form parser with automatic sorting and errorchecking. Tabular input has several advantages for the user, includingdata transparency, the option to maintain standard material libraries,efficient input of digitized or published data and easy electronicexchange of models.

Fixed values or tables may be applied to the following quantities inEquations 1 and 4: ∈_(r), σ, C_(p), k and W_(b). Physical quantities areassigned to materials. The program accepts up to 127 different materialtypes. The elements of the mesh are divided into regions that representphysical objects like tissue types, dielectrics or electrodes. In turn,an element with a given region number is assigned an initial materialand may be associated with different materials, depending on its thermalhistory. Flexible material assignment is the key to the implementationof phase changes. The properties of a region represented by a singletable are non-linear but reversible. If the tissue cools, it returns toits original state.

The algorithm according to present disclosure represents permanent statechanges by reassigning the material reference when the element reaches asetpoint temperature. The algorithm also maintains values of theArrhenius damage integral that represent the extent of chemical changes.Damage integral values can also be used as setpoints. With thereassignment method, properties can be permanently changed, even if theelement cools. For example, a particular element 14 is originallyassociated with a material with electrical conductivity that dropssharply at 100° C. At this temperature, the reference is shifted to amaterial with the properties of desiccated tissue. The conductivity willnot rise even if the electrical power and temperature drops. Allelements of a region have the same region number but may have differentmaterial associations. A region may have up to five alternate materialswith ascending temperature or damage-integral setpoints.

5. Fixed-Current Operation and Temperature Control

As mentioned in Section 3, the RF solution yields the quasi-staticpotential at all points in solution volume for specified electrodepotentials. In many applications, the equipment operates at fixedcurrent or power rather than voltage. These global quantities are notnatural boundary conditions of Equation 1, the method of the presentdisclosure uses an indirect approach. The current or power is calculatedusing the present voltages and material states. An adjustment factor toachieve the desired state is calculated and used to scale potentialvalues.

The global power P₀ is calculated by taking a volume integral oftime-averaged resistive power density (σE*/2) over all material elementsof the solution volume. In multi-electrode systems, the method assigns ageneral amplitude change to all electrodes. If P is the desired powerlevel, then potential values are scaled by √{square root over (P/P₀)}.The calculation of current from an electrode requires a surface integralof current density over the fixed-potential region. Initially, themethod gathers a list of facets of fixed-potential region elements thatare common to conductive/dielectric elements. The program determines thecurrent passing through each facet using a 4×4 normal coordinateintegral of normal electric field, and sums the real and displacementcurrents to find the amplitude I₀. If I is the desired currentamplitude, potential values at all nodes are scaled by I/I₀.

The tabular functions described in Section 4 are also employed in themethod of the present disclosure to represent arbitrary temporalvariations such as a time-dependent thermal source q_(m) or boundarytemperature. The program also stores several generalized functions thatare used for feedback control. In this case, the method uses aregion-averaged temperature as the independent variable and use thetabular function output as a multiplying factor for probe voltages. Themethod generates a programmed voltage to maintain average temperature ina region near 60.0° C. The control function capability makes the methodof the present disclosure a useful tool to test the dynamics andstability of control loops in electrosurgery equipment.

6. Sequence Control

Because of wide variety of operating modes encountered in electrosurgeryand RF ablation procedures, the system and method of the presentdisclosure include a versatile control capability based on a custominterpreter language. The program makes decisions and change operatingmodes depending on the temperature history of the medium. Theorganization of the control script follows the logic of an actualtreatment. Initially there are given conditions that are not controlledby the operator: the geometry of the patient and probe and the intrinsicproperties of materials (including temperature dependencies). During theprocedure, the operator modifies applied voltages and the pulse sequencein response to changes of electrical and thermal properties.Accordingly, the control script contains one SETUP section for thedefinition of general control parameters and material properties. Thissection is followed by an unlimited set of STATE sections. Each STATEsection represents a particular activity that the operator may choose,such as, for example, “apply 2000 mA current until the probe impedancerises 10% from its initial value,” or “lower the voltage to half itspresent value until the maximum temperature of the region drops to lessthan 60° C.”

STATE commands perform the following exemplary types of activities: 1)setting STATE properties such as the maximum duration and frequency ofcommand updates; 2) setting parameters of the RF solution such asvoltage amplitude and phase, frequency or global power; 3) returningvalues of state variables such as average temperature in a region,elapsed time in a STATE, and electrode current; 4) storing numbers inprogram variables with the operation to perform a broad range ofmathematical operations; 6) testing conditions to exit a STATE; and 7)implement control loops to repeat a group of STATES.

Table 2 below shows an example of the control language.

TABLE 2 X(4) = 28.0 Relax X(10) = Impedance(3) X(2) = X(10) 1.15 *ENDSTATE #STARTLOOP STATE 200.0 2.5 Voltage(3) = X(4) Relax ShowImpedance(3) Show MaxTemp X(3) = StateTime X(1) = Impedance(3) TestX(1) > X(2) ENDSTATE STATE 0.0 0.0 Test X(3) > 10.0 X(5) = X(4) 1.0 -X(4) = X(5) ENDSTATE STATE 15.0 15.0 Voltage(3) = 2.0 Relax ENDSTATE#GOTO STARTLOOP 402.0

The algorithm according to the present disclosure records a variety ofdiagnostic data during the thermal solution. The program creates datafiles (snapshots of the RF or thermal solutions) at up to 500 specifiedtimes. The program also allows the user to place up to 20 diagnosticprobes at different positions in the solution volume. The probesgenerate output files containing complete records of temporal variationsof electric and thermal quantities.

7. Analysis Tools

The method of the present disclosure includes dedicated post-processorsto generate plots and to analyze RF and thermal solutions. The methodfeatures interactive graphical capabilities that allow the user to moveeasily through the solution volume. Two dimensional plots showquantities in a plane normal to the x, y or z axes. There are twostyles. First, is a plane plot which is a projection of field quantitiesto a rectangular grid. The simple grid makes it possible to constructadvanced displays such as 2D and 3D filled-contour plots. Second, is aslice plot which is a projection that preserves the structure of theconformal mesh. This mode accurately shows regions boundaries and allowspoint-and-click display of point or scan quantities in a particularplane of the medium 8.

An exemplary slice plot 30 of the medium 8 is shown FIG. 5. The sliceplot 30 displays an electrosurgical probe 32 generating an electricfield 33 represented by contour lines 34. RF solution quantities thatare displayed include the quasi-static potential, electric fieldcomponents, real and displacement currents and power density. Theprogram shows either the amplitude of field quantities or a snapshot ofthe field at a given reference phase.

The method according to the present disclosure also plots spatialvariations of ∈_(r)(T) and σ(T) for solutions with temperature-dependentmaterials. Thermal quantities plotted include the temperature andthermal flux. Spatial variations in k(T), C_(p)(T), and W_(b)(T) aredisplayed for solutions with temperature-dependent materials.

Solutions are also represented using advanced 3D plots with surfacesconstructed directly from the conformal mesh. FIGS. 6A-B show twoexamples. FIG. 6A shows an isothermal surface with identifying regionsthermal flux of the medium 8. In particular, FIG. 6A shows the probe 32having an exposed electrically conductive shaft 31 and an insulatedsheath 35 within a tissue 36. It is envisioned that color coding (shownin gray scale) may be used to illustrate gradual shifting of the thermalsolution. On the right side of the plot, a gray scale coded index 40 isshown. It is envisioned that the index 40 may also be color-coded tobetter represent the gradations in thermal flux of the plot. This typeof plot is useful for displaying the shape of the treatment volume.

FIG. 6B shows boundary facets on selected regions. FIG. 6B also showsthe probe 32 having with the shaft 31 and the sheath 35. This plot typeis useful to display surface temperatures on electrodes and organs. Bysplitting objects into two or more regions, it is possible to displaytemperature variations on arbitrary three-dimensional surfaces. On theright side of the plot, a gray scale coded index 42 is shown. It isenvisioned that the index 42 may also be color-coded to better representthe gradations in temperature of the plot.

Both post-processors have extensive capabilities for quantitativeanalysis. These include point calculations and line scans using asecond-order, least-squares-fit procedure. Information is displayedwithin the program or is ported to an analysis history file in textformat. Line scans generate screen plots with digital oscilloscopefeatures. Complex or repetitive operations are controlled by scripts.The program performs automatic volume integrals over regions of thesolution space. Calculated quantities include RF field energy, totalresistive power dissipation, and average temperature. Surface integralsemploy the same routines for current calculations. Output quantitiesinclude net current or thermal flux. The post-processors include theoption to create matrix files where electrical or thermal quantities areprojected from the conformal mesh to regular box mesh.

8. Application Examples

This section reviews two examples that illustrate some of the featuresof the system and method of the present disclosure. The first exampleaddresses solution accuracy, operation at constant current and theability to handle temperature-dependent materials. The simple testgeometry consists of concentric cooled spherical electrodes of radiir_(i)=1.0 cm and r_(o)=5.0 cm separated by a uniform medium with∈_(r)=1.0. The electrical conductivity is defined by a table with valuesthat start at σ=1.0 S/m at T=0.0° C. and drop smoothly to 0.2 S/m in therange T=20.0° C. to T=120.0° C. Both electrodes are maintained at areference temperature of 0.0° C. The outer electrode is grounded and theinner electrode has a time-dependent voltage to maintain constantcurrent. The medium has the thermal properties k=100.0 W/m-° C.,C_(p)=4800.0 J/kg-° C. and ρ=1000.0 kg/m³. The approximate side lengthof the elements near the inner sphere is 0.1 cm.

The total current flowing through the medium is determined by a surfaceintegral of normal electric field over the inner sphere. Initially allelements are near 0.0° C. and the conductivity is σ≈1.0 S/m. This allowsfor comparison of the impedance calculated by the algorithm to theanalytic formula:

$\begin{matrix}{R = \frac{\left( {\frac{1}{r_{1}} - \frac{1}{r_{0}}} \right)}{4\; \pi \; \sigma}} & (5)\end{matrix}$

Inserting values in Equation 5 gives 6.366Ω, while the calculation usingthe method according to the present disclosure yields 6.371Ω. The smalldifference results from the fact that the shape used by the method ofthe present disclosure is a multi-faceted polyhedron rather than asphere, therefore it is a close approximation of the actual shape of themedium 8.

The thermal solution advances for 50 seconds with a variable time step.The control script contains three STATES. The first sets up a baselineRF solution at f=60.0 Hz with 250 V applied to inner sphere. The secondSTATE advances the thermal solution for 10 seconds and updates commandsat 1 second intervals to include effects of changing electricalconductivity. The STATE includes the command Current(3)=5.7.

The command maintains approximately constant current from the innersphere by correcting the voltage. In response to the command, thealgorithm recalculates the field every 1.0 second, takes the surfaceintegral to determine the current, and scales potential values in thesolution space to bring the current to 5.7 Amps. Additional commandsinstruct the program to record values of voltage and maximum temperaturein the solution volume at each update. The third state (which runs from10 to 50 seconds) has the same commands with an update interval of 2.5seconds.

FIG. 7 shows the voltage variation required to maintain constant currentwith decreasing conductivity for two values of the drive current. The5.8 Amp curve 50 exhibits thermal run-away, while a curve 52 does nothave a similar defect, A constant-current solution may become unstablewhen the electrical conductivity drops with increasing temperature.Thermal run-away is a real effect rather than a numerical artifact;similar phenomena are observed in experiments with electrosurgeryequipment. FIG. 7 demonstrates that non-linear solutions must beapproached with caution—a small change in input may make a dramaticdifference in output.

FIG. 8 shows a two dimensional slice plot of the electrical solution.Here, the probe 32 with length 3.0 cm is inserted in a 2.0 cm diametertumor 54 within a liver 56. The calculation was performed in singleoctant with symmetry boundaries at x=0.0, y=0.0 and z=0.0 using a meshof 74,000 elements. In the calculation, both the liver 56 and tumor 54have the conductivity variation shown in FIG. 4. Electrical conductivityrises 2% per degree up to around 100° C. where a sharp drop representstissue vaporization.

Table 2 above shows the main portions of the STATE sections of thecontrol script. The included commands from the initial STATE store avalue of applied voltage in program variable X(4) as a parameter,calculate a baseline normalized solution to determine thezero-temperature impedance and set a trigger level 15% higher. The nextthree STATES are contained in a loop that runs for a maximum time of 402seconds. The commands of the first STATE sets the voltage to X(4) andadvance the thermal solution until either: 1) the state time reaches 200seconds or 2) the impedance exceeds the trigger level. The RF solutionis updated and checks are performed at 2.5 second intervals.

The second STATE has zero duration so the commands are executed onlyonce per loop cycle. No action is taken if the duration of the previousstate exceeds 10 seconds. Otherwise, the baseline voltage is reduced by1.0 V. The commands of the third STATE reduce the voltage to 2.0 V andadvance the thermal solution for 15 seconds to allow cooling of theregion near the probe.

The algorithm according to the present disclosure also outputs data in avariety of graph, such as impedance over time and temperature over timeplots shown in FIGS. 9A and 9B respectively. In the first pulse cycle,the impedance initially drops as the medium heats and then follows arapid rise as the medium around the probe tip approaches 100° C. In theimpedance plot, triggering does not occur exactly at the setpoint levelof 300° because checks are performed only at 2.5 second intervals. Thetemperature plot of FIG. 9B shows the maximum temperature recorded inthe solution space at a point near the tip of the probe 32 in the regionof enhanced electric field. Records of temperature at other positionsconfirm that the pulse cycle restricts temperature in the bulk of mediumsurrounding the probe to values less than 100° C.

In addition, the algorithm also displays electrical conductivity plots.FIG. 10 shows the spatial distribution of electrical conductivity in theplane y=0.0 late in the pulse (t=65.1 seconds). Although some areas haveenhanced conductivity, σ is sharply reduced near a region 60representing the tip of the probe 32. On the right side of the plot, agray scale coded index 62 is shown. It is envisioned that the index 62may also be color-coded to better represent the gradations in electricalconductivity of the plot.

The simulation involves 120 recalculations of the electrical field witha total run time of 22 minutes on a 3 GHz computer. The exampleillustrates that complex three-dimensional simulations are feasible onstandard personal computers with careful pre-analysis to eliminateunnecessary details.

FIG. 11 shows a flow chart illustrating the method according to thepresent disclosure. In step 300, the processor 102 generates the medium8. The medium 8 may include the tissue 20 with the probe 32 therein asshown in FIG. 3. In step 302, the processor 102 divides the medium intohexahedron elements 14 to generate the three-dimensional conformal mesh12 as shown in FIG. 2C. In step 304, the processor 102 assigns eachelement 14 to a specific material region. For example, elements 14within the protrusion of the probe 32 as shown in FIG. 3 are assigned toa material region of the probe, while elements 14 within the tissue 20are assigned into the corresponding material region. In step 306, theprocessor 102 calculates the electrical energy solution and the thermaltransport solution. In step 308, the display device 112 outputsgraphical representations of the solutions, such as various plotsillustrated in FIGS. 5-10.

While the disclosure has been shown and described with reference tocertain preferred embodiments thereof it will be understood by thoseskilled in the art that various changes in form and detail may be madetherein without departing from the spirit and scope of the disclosure.

1. A method for performing a simulation of an electrosurgical procedure,the method comprising the steps of: generating a three-dimensional modelof a medium including at least one material region representing at leastone of a physical object and a region of the tissue; calculating atleast one of an electrical energy solution and a thermal transportsolution based on at least one mathematical model and finite elementmethod; and updating at least one property value of the at least onematerial region at least once in response to temperature changes basedon an empirically derived value, wherein at least one of the electricalenergy solution and the thermal transport solution is updated as afunction of the at least one property value.
 2. The method according toclaim 2, further comprising the steps of: dividing the model of themedium into a plurality of hexahedron elements to generate athree-dimensional conformal mesh; assigning each element to acorresponding at least one material region; and calculating at least oneof an electrical energy solution and a thermal transport solution basedon at least one mathematical model, the plurality of hexahedron elementsand finite element method.
 3. The method according to claim 2, furthercomprising the step of: generating and displaying at least one graphicalrepresentation representing at least one form of the electrical energysolution and one form of the thermal transport solution on a outputdevice.
 4. The method according to claim 1, wherein the step ofcalculating the electrical energy solution further comprises a step ofsolving equations 1 and 2: $\begin{matrix}{{\nabla{\cdot \left\lbrack {\left( {ɛ_{r} - \frac{j\sigma}{ɛ_{0}\omega}} \right){\nabla\varphi}} \right\rbrack}} = 0} & (1) \\{E = {- {\nabla\varphi}}} & (2)\end{matrix}$ wherein ∈_(r) is a relative dielectric constant, and σ iselectrical conductivity of each of the elements, j is an imaginarynumber, ∈₀ is vacuum permittivity, ω is angular frequency, φ is aquasi-static potential, which is a complex number, and E is amplitudeand phase of electric field.
 5. The method according to claim 4, whereinthe step of calculating the thermal transport solution further comprisesa step of solving equations 3 and 4: $\begin{matrix}{q_{e} = {\sigma \; {E \cdot E}*{/2}}} & (3) \\{{{pC}_{p}\frac{\partial T}{\partial t}} = {{\nabla{\cdot \left( {k{\nabla T}} \right)}} + q_{e} + q_{m} + {W_{b}\left( {T_{b} - T} \right)}}} & (4)\end{matrix}$ wherein ρ is density of the medium, Cρ is specific heat, kis thermal conductivity, qm is heat generated by metabolism, qe isresistive source, and Wb(Tb−T) represents thermal variations caused byblood flow.
 6. The method according to claim 1, wherein the at least oneproperty value is representative of at least one physical propertyselected from the group consisting of relative dielectric constant,electrical conductivity, specific heat, thermal conductivity and bloodflow.
 7. The method according to claim 6, wherein the at least onephysical property is temperature-dependent.
 8. The method according toclaim 6, wherein a tabular function stores a plurality oftemperature-dependent values for the at least one physical property. 9.The method according to claim 2, wherein the plurality of elements areordered in the three-dimensional mesh according to cubic logic.
 10. Themethod according to claim 2, wherein the graphical representation is athree-dimensional plot and the at least one form of the thermal solutionis an isothermal surface of the medium.
 11. The method according toclaim 2, wherein the graphical representation is a two dimensional sliceplot and the at least one form of the electrical solution is electricalconductivity of the medium.
 12. A system for performing a simulation ofan electrosurgical procedure on tissue, the system comprising: at leastone processor coupled to a computer-readable media configured to store aset of instructions capable of being executed by the at least oneprocessor, the instructions including the steps of: generating athree-dimensional model of a medium including at least one materialregion representing at least one of a physical object and a region ofthe tissue; calculating at least one of an electrical energy solutionand a thermal transport solution based on at least one mathematicalmodel and finite element method; and updating at least one propertyvalue of the at least one material region at least once in response totemperature changes based on an empirically derived value, wherein atleast one of the electrical energy solution and the thermal transportsolution is updated as a function of the at least one property value.13. The system according to claim 12, wherein the instructions furthercomprise the steps of: dividing the model of the medium into a pluralityof hexahedron elements to generate a three-dimensional conformal mesh;assigning each element to a corresponding at least one material region;and calculating at least one of an electrical energy solution and athermal transport solution based on at least one mathematical model, theplurality of hexahedron elements and finite element method.
 14. Thesystem according to claim 12, further comprising: an output deviceconfigured to output at least one graphical representation representingat least one form of the electrical energy solution and one form of thethermal transport solution.
 15. The system according to claim 12,wherein the at least one property value is representative of at leastone physical property selected from the group consisting of relativedielectric constant, electrical conductivity, specific heat, thermalconductivity and blood flow.
 16. The system according to claim 15,wherein the at least one physical property is temperature-dependent. 17.The system according to claim 15 wherein a tabular function stores aplurality of temperature-dependent values for the at least one physicalproperty.
 18. A computer-readable storage medium storing a set ofcomputer-executable instructions for performing a simulation of anelectrosurgical procedure, the computer-executable instructionscomprising the steps of: generating a three-dimensional model of amedium including at least one material region representing at least oneof a physical object and a region of the tissue; calculating at leastone of an electrical energy solution and a thermal transport solutionbased on at least one mathematical model and finite element method; andupdating at least one property value of the at least one material regionat least once in response to temperature changes based on an empiricallyderived value, wherein at least one of the electrical energy solutionand the thermal transport solution is updated as a function of the atleast one property value.
 19. The computer-readable storage mediumaccording to claim 18, wherein the computer-executable instructionsfurther include the steps of: dividing the model of the medium into aplurality of hexahedron elements to generate a three-dimensionalconformal mesh; assigning each element to a corresponding at least onematerial region; and calculating at least one of an electrical energysolution and a thermal transport solution based on at least onemathematical model, the plurality of hexahedron elements and finiteelement method.
 20. The computer-readable storage medium according toclaim 18, wherein the computer-executable instructions further includethe step of: generating and displaying at least one graphicalrepresentation representing at least one form of the electrical energysolution and one form of the thermal transport solution on a outputdevice.